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INTRODUCT ION 


This  technical  report,  which  is  the  last  in  a series,  covers  new 
results  in  linear  restoration  of  blurred  photon-limited  images.  Our  most 
recent  previous  technical  report  [1]  is  highly  recommended  as  background 
for  this  report,  and  we  shall  freely  call  upon  results  derived  there. 
However,  this  report  does  contain  a brief  review  (together  with  some 
minor  corrections)  of  the  earlier  work. 

For  the  purposes  of  motivation,  we  should  mention  at  the  start  why 
wo  are  interested  in  space-variant  linear  restoration  of  degraded  images. 
Compensated  imaging  systems  correct  for  atmospheric  blur  over  only  a fi- 
nite field  (the  so-called  "isoplanatic  region").  If  objects  larger  than 
this  field  are  encountered,  the  blur  present  in  the  image  will  be  markedly 
space-variant.  Furthermore,  while  conventional  linear  restoration  theor- 
ies assume  stationary  object  statistics,  in  practice,  real  ensembles  of 
objects  are  usually  nonstationary  in  one  or  more  attributes.  Both  space- 
variant  blurs  and  nonstationary  objects  l^ad  to  space-variant  filtering 
as  a desired  approach. 

Section  II  to  follow  contains  the  review  mentioned  above.  In  Section 
III,  we  describe  an  algorithm  for  finding  a linear  restoration  filter 
that  maximizes  the  particular  image  quality  criterion  adopted  here.  Sec- 
tion IV  gives  the  results  of  a study  of  restoring  photon-limited  images 
with  linear  shift-variant  restoration  filters. 


II . EQUATIONS  TOR  THE  DISCRETE  SPACE-VARIANT  RESTORATION  PROBLEM 

As  stated  above,  in  this  section  we  give  a review  of  the  discrete 
notation  and  equations  used  in  our  study  of  space-variant  restoration 
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of  photon-1 imi tod  images.  During  this  review,  we  will  show  some  correc- 


tions to  the  equations  given  in  Section  IV  of  reference  1. 


(a)  Discrete  Notation 

The  physical  process  of  the  formation  of  an  image  is  modeled  by  a 

linear  matrix  expression.  The  continuous  object  radiance  distribution  is 

represented  by  a column  vector  o whose  elements  are  M equally  spaced 

- o 

samples.  The  use  of  lexicographic  ordering  to  represent  two-dimensional 
object  and  image  distributions  is  discussed  in  reference  1. 

The  formation  of  an  image  from  the  object  is  modeled  by  the  equation 


i = [B] 


(1) 


Here,  the  image  intensity  distribution  is  represented  by  the  elements 
of  the  column  vector  _i.  The  blur  matrix  [B]  consists  of  M samples 
of  each  of  the  impulse  responses  of  the  imaging -blurring  system.  In 

general,  the  number  of  samples  of  the  image  might  not  equal  the  num- 
ber of  samples  of  the  object  M . 

o 

The  detected  image  vector  d consists  of,  again  in  general,  M ele- 

~ d 

ments , each  of  which  represents  the  photocount  from  one  element  of  a dis- 
crete array  of  photodetectors.  In  our  work,  we  have  assumed  that  M = M 

d i 

and  that  the  n**1  element  d^  of  d de-pends  directly  on  only  the  n^*1 

element  i^  of  _i.  In  fact,  we  model  the  image  detection  process  by 

assuming  that  d is  Poisson  distributed  with  mean  rate 
n 


X -^i  , 

n h'v  n 


(2) 
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where  A is  the  area  of  a detector  element,  and  T is  the  observation 
interval . 

The  detected  image  d is  restored  by  linear  filtering  with  the 
restoration  filter  [H].  This  filter  will  be  chosen  to  make  the  restored 
image  r,  given  by 

r = [H]  d , (3) 

as  "close"  as  possible  to  an  "ideally  filtered"  object 

o = [S]  o . (4) 

Here,  [s]  is  an  M~  X Mq  matrix  of  samples  of  the  impulse  response  of 
an  ideal  filter.  The  normalization  indicated  by  A was  defined  in  refer- 
ence 1.  In  general,  therefore,  the  filter  [H]  is  of  dimension  M~xM  . 

o d 

(b)  Minimum  Mean-Squared-Error  Filter 

One  method  of  making  the  restored  image  r "close"  to  the  ideally 
filtered  image  o is  to  minimize  the  mean-squared  value  of  the  error 
column  vector  e given  by 

e = o - r = [s]  o - [H]  d . (5) 

Thus,  we  wish  to  minimize  the  quantity 

6 = E^e  | = EjTrUeSj  (6) 

where  t signifies  matrix  transpose,  E is  the  statistical  expectation 
operator,  and  Tr(  ) is  the  matrix  trace  operation. 

Referring  to  reference  1 for  the  derivation  of  the  following  re- 
sults, the  meun-squured-error , in  general,  is  given  by 
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An  alternate  form  of  the  mean-squared-error  actually  used  in  calculations 
is 


Tr<k2N2[S]  [Sq]  [S]*  +N|[H]| 


' A _ 

[B]  o 


+ N[B][K  ] [B ] 1 J - 2kN[S][K  ] CB]t  ] [H]1" 
o o 


(8) 


The  appropriate  minimum  mean-squared -error  filter  can  be  found  in  the 
literature  [2]  and  is  given  by 


^MMSE1 


kN[S][<R  ] [B]1 
o 


/r 

\L« 


S A.  2 

[B]  o 


^-1 


+ N[B][K  ] [B]1" 


(9) 


This  filter  will,  in  general,  be  space  variant  and  nonsquare,  depending 

A A ^ V 

on  the  properties  and  sizes  of  [S],  [B] , o,  and  [JR  ] . 

In  order  to  study  the  performance  of  a given  restoration  filter, 

some  measure  of  imago  quality  must  be  defined.  ’.Ve  have  chosen  to  define 
image  quality  Q by 
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Q = 


where  6 is  the  mean-squared-error  and  & is  the  expected  signal  en- 
ergy at  the  output  given  by 


A N2Tr  (l 


] [B 1 [H  ] 
o 


*)  • 


Substitution  of  the  minimum  mean-squared -error  filter  [H  ] into  one 

ARISE 

of  the  expressions  for  mean-squared-error,  Eq.  (7)  or  (8),  yields  the 
expression  for  the  minimum  possible  mean-squared-error  with  linear  res- 
toration 


6 . = k2N2  Tr  [S][fR  J [i] 

mm  o 


[B]t  LB]  5 


+ N[B]  [£R  ] [B]^  [BHtf]  CS^ 
o o 


(In  reference  1,  the  corresponding  equation,  Eq.  (36),  has  a misplaced 
parenthesis.)  In  our  computations,  the  alternative  form 


2-2  / a v * t 1 

6 . = k N Tr  (fSHa  Us]  - i 

min  \ o k 


BLn-JCBlfe  1 CS]1 
MMSE  o 


for  the  minimum  mean-squared-error  has  been  used.  (The  corresponding 
Eq.  (37)  in  reference  1 has  an  omission.  Note  here  that  the  form  of  Eq. 
(38)  in  reference  1 used  for  comparison  to  the  continuous  case  is  cor- 


rect only  if  all  matrices  in  the  expression  are  square,  i.e.,  if  1Uq  = 

11.  = M~.  The  more  generally  valid  expression,  requiring  only  M = 
i o i o 
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2—2  * ^ t 

5 . = k N Tr  [S]  [fR  ] [S] 

mm  o 


[Bio  +K[B][«][B]t 

— o 


'[B]  5 


which  in  the  square  matrix  case  reduces  to  Eq.  (38),  reference  1.) 

Thus,  we  have  several  forms  for  the  denominator  of  the  expression 
for  Q.  For  the  numerator,  we  substitute  Eq.  (9)  into  Eq . (11),  yielding 


W = k 5 Tr  tsn^lCB]* 


[B]  o + N [B]  [£R  ] [B]t 

— o 


A V /\  f- 

• [B]  [9?  ][B]r 


* * - - A ^ -f-  \ A V/  /\  + \ 

[B]  o + N CB ] C£Rq  ] [B]  [BHSUCsr  . 


(In  reference  1,  Eq.  (43)  contains  an  exponent^ error  and  Eq.  (44)  is  more 
correctly  written 


MMSE  S 


for  the  minimum  mean-squared-error  Q.) 

As  we  show  in  Section  III,  we  can  improve  on  the  performance  pre- 
dicted by  Eq.  (16).  Specifically,  we  can  do  better  than  the  minimum 
mean-squared— error  filter  by  using  a filter  which  maximizes  the  quality 
factor  Q. 
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III.  maximum  image  quality  kilter 

Equation  (1C)  is  not  the  maximum  obtainable  Q using  linear  res- 
toration techniques.  This  is  true  because  in  minimizing  the  mean-squared- 
error  we  are  not  necessarily  maximizing  the  ratio  of  signal  to  mean- 

squared  error.  Tims,  we  would  like  to  find  the  filter  [H  ] which 

MQ 

maximizes  the  image  quality 


Q A i . 

6 


(10) 


By  following  a procedure  similar  to  that  used  for  finding  the  minimum 
mean-squared -error  filter  (see  Appendix  A),  two  necessary  conditions  are 
found  for  any  fH^  ] producing  a local  maximum  in  Q: 


if  Tr  (2  [Hmq]  [B]  BU  Cb]*  [G]t)  Tr  (2^  ] [D]  fc] 


2kij2[G][B][K  Hs]*) 

O ! 


(17) 


and 


2NQ  Tr  [G] 
max  ' 1 


[B]  o 


+ 5 to]  [a  1 Cb]1 

Vx 


CG] 


< 0 


(18) 


for  all  matrices  [G] , where  in  (17)  expressions  (7)  and  (11)  have  been 
used  and 


[D] 


+ N[Bm  Hb]1 

o 


(19) 
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By  rearranging  terms,  we  find  that  condition  (17)  will  be  satisfied 


if  and  only  if 

A \ 1 

[B] [£*  1 [D]t  (20) 

° / 

and  condition  (IB)  if  and  only  if  (20)  and 

NCBHSR  HbI*  > [0]  (21) 

Qmax 

are  satisfied,  where  [0]  is  a matrix  of  all  zero  elements.  The  ine- 
quality (21),  for  our  purposes,  is  simply  an  added  condition  which  must 
be  satisfied  by  any  satisfying  (20).  It  is  no  help  to  us  in 

finding  our  maximum  Q filter.  TJote  that,  if  the  factor  (Q  -1)/Q 

max  max 

is  regarded  as  a simple  constant,  the  maximum  Q filter  differs  from 
the  minimum  mean-squared -error  filter  only  by  this  simple  internal  pa- 
rameter . 

Equation  (20)  is  not  a true  solution  for  the  elusive  maximum  Q 

filter  since  it  is  necessary  to  know  a priori  the  value  of  Q , and 

max 

in  our  studies  we  need  [H,  1 to  find  Q . In  a real  image  restora- 

MQ  max 

tion  system,  however,  this  might  not  be  a problem,  as  [H^l  might,  for 

example,  be  determined  experimentally  from  the  minimum  mean-squared-error 

filter  by  manually  adjusting  the  factor  (Q  - 1)/Q  . It  should  be 

max  max 

mentioned  here  that,  not  only  is  it  difficult  to  find  a filter  satisfying 
(20)  and  (21),  but  these  arc  only  necessary  conditions  for  the  maximum 
Q filter.  These  conditions  guarantee  only  a local  maximum  in  the  Q-fH] 
space . 
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[Hi  = kNlsns  ItB]* 

MQ  o 


[Bl  o 


+ N 


/Q  - 1' 
' max 


ax 


For  tho  purposes  of  this  study  we  have  developed  a procedure  for 


finding  the  locally  maximum  Q filter 


V • 


The  procedure  may  be 


described  as  solving  the  following  pair  of  recursive  equations: 


- [Hi  = [S][£R  ItS]*1 
kN  1 


- /\  £ 
[B]  o 


♦ sfei ')  [bj  ] [e  ] 1 

\ Qi-i  / 


and 


-2 
= N 


Tr(-4  [HJCBU&  ItB]*  i [H  lM 
\kN  1 ° kR  i / 


Tr  |[S]  [9^3  [S] 


* + n(4  [H 

VkN  1 


HD]  - 2 IS]  [S  llBl* 
o 


/ kN  1 y 


(22) 


(23) 


The  solution  is  reached  by  performing  these  steps: 

(1)  Choose  an  initial  value  Q by  some  method,  (For  this 

o 

Q^,  one  might  use  the  value  of  Q obtained  by  using 
the  space-variant  or  space-invariant  minimum  mean- 
squared-error  filter,  or  one  might  try  several  values 
for  and  choose  that  value  which  produces  the  larg- 

est . Our  experience  is  that,  in  most  cases,  a value 
of  Qq  slightly  more  than  unity  works  well.) 

(2)  Find  [H  ] from  Eq.  (22). 

(3)  Use  [Hj]  to  compute  using  Eq.  (23). 

(4)  Repeat  steps  2 and  3 always  using  the  latest  for  the 

next  cycle. 

(5)  Stop  when  the  values  of  are  no  longer  increasing. 

(6)  Check  to  insure  a maximum  rather  than  an  inflection  point 

by  entering  Eqs . (22)  and  (23)  with  several  values  of 

<3  on  both  sides  of  the  final  Q.  from  step  4. 
i-1  l 
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Wo  have  :‘ound , using  this  procedure,  solutions  for  the  maximum  Q 
filter  much  more  easily  than  we  had  expected  when  the  procedure  was  first 
proposed.  Also,  although  we  have  no  guarantee  that  we  have  found  the 
global  maximum  in  the  Q-Dl]  space,  we  have  tested  a wide  range  of  val- 
ues of  the  factor  (Qq-1)/Qo  and  have  never  found  a Q higher  than  that 
obtained  from  the  maximum  Q procedure. 

IV.  LINEAR  SPACE-VARIANT  FILTERING  RESULTS 

In  this  section,  we  present  results  of  computer  analysis  of  restored 
image  quality  possible  when  using  linear  space-variant  restoration  tech- 
niques. We  have  used  both  minimum-mean-squared-error  and  maximum-Q-algo- 
rithm  filters.  We  begin  by  describing  the  notation  and  some  of  the  object 
statistics  and  blur  types  used  through  Section  iv. 

(a ) Computation  Parameters 

In  our  computer  studies,  we  have  treated  several  types  of  nonsta- 
tionary objects  and  space-variant  blurs.  We  usually  allow  either  the 
object  statistics  or  the  blur  impulse  response  to  be  position  dependent, 
but  not  both  at  once.  Later,  we  will  define  the  exact  forms  used  for 
those  statistics  and  impulse  responses.  At  this  point,  we  shall  describe 
some  general  assumptions  and  define  the  types  of  space-invariant  blur 
and  stationary  object  used  throughout. 

First,  in  all  cases,  we  have  restricted  ourselves  to  one-dimensional 
objects.  During  the  period  covered  by  this  contract,  the  speed  and  stor- 
age capabilities  of  the  computer  at  our  disposal  were  never  sufficiently 
large  to  allow  two-dimensional  image  processing.  In  fact,  in  nearly  all 
of  our  studies,  we  restricted  the  image  and  object  sizes  to  35  pixels. 
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In  only  a few  cases  did  we  allow  object  and  image  sizes  to  grow  to  49 
pixels.  In  these  cases,  we  were  checking  for  object-dependent  changes 
in  image  quality  due  to  edge  effect  errors.  As  the  object  size  grew, 
all  the  Q curves  shifted  upward  slightly  but  always  by  nearly  the  same 
amount . 

As  we  stated  in  Section  II,  the  sizes  of  the  object,  classical  in- 
tensity image,  detected  image,  and  ideal  image,  M , M , M, , and  M~, 

o i d o 

respectively,  will,  in  general,  be  different.  This  might  be  true  due  to 
the  physical  realities  of  the  imaging  system  or  because  it  might  be  ad- 
vantageous to  choose  them  so.  If,  for  instance,  the  impulse  responses 
of  the  actual  blur  and  ideal  blur  are  assumed  to  be  spatially  limited, 
the  sizes  of  the  vectors  and  matrices  are  often  chosen  to  retain  all 
available  information  without  suffering  any  unneeded  storage  penalty. 

In  our  studies,  we  have  assumed  more  "realistic"  blurs  which  are 
band  limited  rather  than  space  limited.  These  must  necessarily  be  arti- 
ficially truncated  at  a rather  arbitrary  point,  especially  with  our 
rather  severe  space  limitations.  We  therefore  have  assumed  that  all 
vectors  have  the  same  dimension  and  that  all  matrices  are  square,  that 
is,  that 

M = M.  = M = = M (24) 

o i d o 


where  M is  the  assumed  matrix  and  vector  dimension.  As  stated  above, 
for  all  cases  presented  here,  M = 35. 

Briefly,  with  respect  to  notation,  if  vector  a has  elements  a^, 
then  the  index  i is  assumed  to  run  from  1 to  M.  If  matrix  [A]  has 


elements  CA  3 . . , 

i j 


then  the  index  i runs  along  a column  and  j runs 


along  a row  from  1 to  M. 
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The  idea1  blur  matrix  used  in  all  our  calculations  is  a sine 

2 

function  sampled  at  the  Nyqulst  rate.  This  sine  impulse  response 
corresponds  to  a diffraction-limited  one-dimensional  aperture  function. 
Explicitly,  the  ideal  filter  matrix  elements  are 


CSJ 


ij 


/sin  («(i  - j)/2)\2 

V «(i  - j>/2  ) ' 


This  equation  is  more  simply  written 


i = j 


[S] 


U " 


U(i  - j)l 


i - j even  and  i ^ j . 
i - j odd 


(25) 


(26) 


In  all  our  calculations,  the  blur  used  is  a Gaussian  blur  with  a 
possibly  posit ion -dependent  width.  Explicitly,  it  is 


where  fxl  means  the  smallest  integer  >x  and  w(x)  = W , a constant, 
whenever  a space-invariant  blur  is  specified.  The  form  of  w(x)  for  a 
space-variant  blur  will  be  described  later. 

In  all  our  work  here,  for  both  stationary  and  nonstationary  object 
ensembles,  we  assume  the  covariance  matrix  is  diagonal.  That  is,  we  as- 
sume the  object  variation  is  uncorrelated  from  point-to-point  in  the  ob- 
ject. For  our  stationary  object  ensembles,  the  autocorrelation  matrix 
is 
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(28) 


[nt  ] = a2 [I]  + k2N2  ll1 
o o — 


2 

where  oq  is  the  variance  of  the  object,  [I]  is  the  identity  matrix, 
and  1 is  a vector  all  of  whose  elements  are  unity.  We  have  two  degrees 
of  freedom  here,  the  variance  and  the  flux  level,  and  In  practice  we 
specify  N and  the  object  squared-mean  to  variance  ratio  A at  the 
center,  defined  by 


(29) 


where  is  the  i**1  element  of  the  mean  object  vector  and  [0q1  is  the 

object  covariance  matrix. 

The  constant  k is  defined  by 


A 

= boqTA 


(30) 


where  b is  the  net  energy  transmittance  of  the  imaging  process.  This 
o 

constant  k has  not  been  specified  explicitly  in  our  calculations  for 
two  reasons.  First,  in  all  our  results,  we  present  image  quality  plotted 
versus  the  mean  number  of  photoevents  per  image  rather  than  versus  the 
total  mean  object  brightness.  Thus,  we  have  no  need  to  translate  total 
brightness  into  total  photoevents  by  using  k.  Second,  we  have  implic- 
itly assumed  that  all  quantities  in  Eq.  (30)  are  deterministic  and  known 
a priori.  Thus,  the  restoration  filter  can  compensate  for  the  effects 
of  the  constant  k and  effectively  remove  k from  our  results. 


13 


(b)  Space -Variant  Blur 


Here,  we  show  results  of  restored  image  quality  for  situations  hav- 
ing stationary  object  ensembles  and  space-variant  blurs.  We  use  minimum 
mean-squared -error  and  maximum  Q algorithm  restoration  filters  of  both 
the  space-invariant  and  space-variant  types. 

The  space-variant  blur  we  use  is  that  of  Eq.  (27)  with 


»(«  = »c 


(31) 


where  W and  W are  the  widths  at  the  center  and  edge,  respectively, 
c e 

A parabolic  width  function  was  also  tested,  and  the  results  were  very 
similar  to  those  presented  here. 

Figure  1 shows  Q vs  log^  ^ f°r  severa^  images  restored  with 
minimum  mean-squared -error  filters.  In  part  '(a),  the  central  object 
mean-to-variance  ratio  A is  1 and,  in  (b) , A is  20.  In  both  parts, 
curves  are  shown  for:  (1)  a space-invariant  blur  with  Wb  = 8 and  with 
its  corresponding  space-invariant  restoration  filter;  (2)  a slightly 
space-variant  blur  with  = 8 and  = 10  with  a space-variant  res- 
toration filter;  and  (3)  a severely  space-variant  blur  of  W = 8 and 

c 

W = 16  also  with  a space-variant  restoration  filter.  These  curves 
e 

represent  the  best  one  can  do  with  linear  minimum  mean-squared-error 
restoration  applied  to  these  particular  blurs.  The  Q values  decrease, 
progressing  from  (1)  through  (3)  in  each  part  because,  for  a given  N, 
the  filter  cannot  restore  an  image  as  well  for  larger  average  blurs. 


1-1 


In  Fig  . 2 , Q is  plotted  for  those  same  blurs  of  Fig.  1,  but  the 
images  have  been  restored  using  a space-invariant  minimum  mean-squared- 
error  filter  with  = 8.  As  in  Fig.  1,  part  (a)  has  A = 1 while 

part  (b)  has  A = 20.  The  most  obvious  feature  of  this  figure  is  the 
greatly  reduced  restored  image  quality  resulting  when  the  space-invari- 
ant filters  are  used  with  space-variant  blurs.  While  the  image  quality 
obtained,  when  the  filter-blur  type  and  width  are  matched  to  the  type 
and  width  of  the  blur,  increases  significantly  with  increasing  N,  the 
image  quality  for  unmatched  filters  actually  approaches  an  asymptote. 

More  generally,  as  can  be  seen  from  Fig.  2 and  other  figures  to 
follow,  we  find  that,  whenever  the  filter  blur  does  not  match  the  actual 
blur,  the  image  quality  is  asymptotic.  For  the  square,  symmetric  blur 
matrices  used  in  our  calculations,  this  asymptote  is  approximately 


Tr  (fc: 1 tV <*i ' (»/ l%t >) _1  &)’  % >) 

Tr^([S][Bf]  - [Sf][B])  ([Bf]2tRof])_1  [B]J  [<Ko][aof] 


(32) 


where  the  subscript  f indicates  matrices  used  to  form  the  particular 

A A 

filter  used.  In  all  our  calculations,  we  choose  [S^]  = [S].  This  will 
always  be  possible  because  the  [S]  used  everywhere  was  our  choice  in 

A A 

the  beginning.  Using  this  choice  of  [S^J  and  assuming  [B^ ] and 

V 

[£R  1 have  inverses  the  asymptote  (32)  becomes 
o 


(33) 


IS 


g.  2.  Q for  restoration  with  a space-invariant  minimum 
mean-squared -error  filter.  Two  A values:  (a)  A = 1 ; 
(b)  A=  20.  Different  blurs:  (1)  space  invariant,  r?  1 
(2)  space  variant,  Wc  = 8 and  W = 10;  (3)  space  vari 


The  most  important  factor  in  this  expression  is  the  "blur  difference" 

A A 

[Bf]  - [Bl.  'if  this  difference  is  large,  the  value  of  the  asymptote 
approaches  unity. 

In  Fig.  3,  Q is  shown  for  a severely  space-variant  blur  of  W£=  8 

and  W =16.  In  part  (a),  A = 1,  and  in  (b),  A = 20.  This  blurred 
e 

image  has  been  restored  using  several  space-variant  and  space-invariant 

filters.  Curves  are  shown  for  minimum  mean-squared-error  filters  matched 

both  to  the  space-variant  blur  and  to  space-invariant  blurs  of  widths 

= 4,  8,  and  12.  Points  indicated  by  small  circles  are  for  restoration 

with  a maximum  Q algorithm  filter.  Crosses  show  restoration  with  a 

space-invariant  maximum  Q algorithm  filter  with  W =8.  There  are 

b 

only  a few  of  these  points  plotted  because  the  algorithm  does  not  guar- 
antee a maximum  unless  it  is  allowed  to  be  space  variant.  Here,  as  in 
most  cases  we  have  encountered,  the-maximum  Q value  falls  approximately 
one  unit  above  the  minimum  mean-squared-error  filter  value  of  Q. 

An  interesting  phenomenon  appears  in  part  (a).  The  space-invariant 
minimum  mean-squared-error  filter  has  performed  slightly  better  than  the 
space-variant  minimum  mean-squared-error  filter  for  small  N.  Later,  in 
Fig.  5,  we  shall  present  an  even  more  pronounced  example  of  this  phenom- 
enon. It  is  not  as  disconcerting  as  it  might  seem  to  see  the  space-in- 
variant filter  perform  better  than  the  space-variant  filter  since,  as  we 
mentioned  in  Section  III,  minimizing  the  mean-squared -error  does  not 
necessarily  maximize  the  image  quality  Q.  In  a sense,  this  "better 
performance"  of  the  space-invariant  filter  is  just  an  accident.  By 
chance,  the  mismatched  blur  has  introduced  an  approximation  to  the  fac- 
tor (Q  - 1)/Q  of  the  maximum-Q-a Igor ithm  filter.  This  mismatch  could 
Just  as  well  have  decreased  the  effectiveness  of  the  filter. 
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g . 3 . Q for  space-variant,  Wc  = 8 and  We  = 16 , blur. 
Two  A values:  (a)  A = 1;  (b)  A = 20 . Several  minimum 

mean -squared -error  filters:  (1)  space  variant,  W£  = 8 
and  We  = 16;  (2)  space  invariant,  Wb  =4;  (3)  space 

invariant,  V/j,  - 8;  (4)  space  invariant,  Wb  = 12.  Two 

maximum-Q-algorithm  filters:  (O)  space  variant,  Wc=8, 
v;e  = 16;  (+)  space  invariant,  v,'b  = 8. 


8 


Figure  4 shows  Q for  a slightly  space-variant  blur  with  V = 

and  W = 10  with  both  space-variant  filtering  and  space-invariant 
e 

W =8  filtering.  Part  (a)  is  for  A = 1,  and  part  (b)  is  for  A=20. 
b 

Circles  are  for  maximum  Q algorithm  filtering.  Note  here  that  the 
space-invariant  filter  results  are  asymptotic  to  a different  level  for 
large  N than  are  the  severely  space-variant  blur  results.  Two  other 
features  come  to  light  if  Fig.  4 is  compared  to  Fig.  3.  First,  as  seen 
in  Figs.  1 and  2,  the  general  performance  here  is  better  than  for  the 
severely  space-variant  blur.  Second,  the  blur  mismatch  here  is  not  suf- 
ficiently large  to  allow  the  space-invariant  minimum  mean -squared -error 
filter  to  produce  better  image  quality  than  the  space-variant  minimum 
mean-squared -error  filter. 

Figure  5 gives  Q for  a space-variant  blur  with  = 16  and 

W =32.  As  before,  part  (a)  has  A = 1,  and  (b)  has  A = 20 . Curves 
e „ 

are  shown  for  a matched  space-variant  filter  and  space-invariant  filters 
matched  to  = 12,  16,  24,  and  32.  Small 'circles  show  filtering  with 
a maximum  Q algorithm  filter.  Note  carefully  that  this  figure  has  a 
different  scale  than  all  others  in  Section  IV.  Here,  the  severity  of 
the  blur  lowers  the  curves  well  below  those  of  the  previous  figures.  As 
mentioned  above,  in  Fig.  5(a)  we  have  a dramatic  example  of  the  space- 
invariant  minimum  mean-squared -error  filters  performing  better  than  the 
space-variant  minimum  mean-squared -error  filter.  Here,  the  severe  mis- 
match of  the  blur  allows  good  accidental  approximation  of  the  maximum  Q 
filter  over  a range  of  N.  In  addition,  the  mismatch  is  enough  to  force 
the  space-invariant  filters’  performances  to  the  large  N asymptote's 
minimum  value  of  unity. 
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Fig.  5.  Q for  space-variant,  VVC  = 16  and  Wc=32,  blur. 
Two  A values:  (a)  A = 1;  (b)  A = 20 . Several  minimum 

mean -squared -error  filters:  (1)  space  variant,  Wc  = 16 
and  Wc  - 32;  (2)  space  invariant,  Wp  = 12;  (3)  space 

invariant,  Wjj  = 16;  (4)  space  invariant,  W)5  = 24;  (5) 

space  invariant,  Wp  = 32.  (G)  Maximum  Q filtering. 
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(c)  Nonstationary  Object 


Here,  we  give  resulting  image  quality  when  nonstationary  object  en- 
sembles are  blurred  with  a space-invariant  blur.  We  study  the  full  range 
of  minimum  mean-squared-error  and  maximum  Q algorithm  filters,  both 
space-variant  and  space-invariant  cases. 

We  consider  two  different  contributions  to  the  nonstationary  sta- 
tistics of  the  object  ensemble:  the  mean  and  the  covariance.  We  allow 
both  of  these  to  be  position  dependent,  individually  and  together.  When 
we  indicate  nonstationary  mean,  we  are  using  o of  the  form 


/i  “ 


oi  = A ( 


W /2 
o 


(34) 


where  W is  an  object  width  parameter  and 
o 


A(x)  = 


1 - x for  |x  < 1 


for  x > 1 


(35) 


For  those  cases  where  we  use  a nonstationary  covariance,  the  form  is 


for  i 4 j 


[0] 


ij 


-fir 

W /2 
o 


(36) 


for  i = j 


In  all  of  Section  IV(c),  the  blur  is  space  invariant  with  W.  = 8. 

b 

Figure  6 shows  Q for  three  different  nonstationary  object  ensem- 
bles with  A = 1.  The  upper  three  curves  are  restored  using  the  appro- 
priate space-variant  minimum  mean-squared -error  filters.  The  lower  three 
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g.  6.  Q vs  logio  N for  different  nonstationary  object  statistics: 
(1)  & (4)  nonstationary  mean;  (2)  & (5)  nonstationary  covariance;  (3)  & 
(6)  nonstationary  mean  and  covariance.  Different  filters  used  are: 
(1)  , (2),&  (3)  space-variant  minimum  mean-squared  error  filters;  (4), 
(5),  & (6)  space-invariant  minimum  mean-squared-error  filters. 


g.  7.  Q for  different  nonstationary  object  statistics:  (1)  nonsta 
tionary  mean;  (2)  nonstationary  covariance;  (3)  nonstationary  mean 
and  covariance. 


curves  correspond  to  restorations  obtained  using  space-invar iant  minimum 
mean-square<|-error  filters.  The  three  curves  in  each  group  include  an 
object  ensemble  with:  (1)  only  the  mean  nonstationary,  (2)  with  only 
the  variance  nonstationary,  and  (3)  with  both  the  mean  and  variance  non- 
stationary. As  can  be  seen,  there  is  little  difference  for  the  three 
types  of  nonstationarity . The  most  pronounced  feature  of  the  curves 
here  is  the  gap  between  the  group  of  three  space-invariant  filter  curves. 
The  nonstationary  mean  curve  is  almost  as  high  as  the  space-variant 
filter  curves.  The  two  curves  showing  space-invariant  filters  used  with 
nonstationary  object  variances  show  definitely  lower  performance. 

Figure  7 has  curves  plotted  for  the  maximum  Q filter  restoring 
the  three  types  of  nonstationarity  listed  for  Fig.  6.  Here,  there  is 
even  less  difference  between  the  performances  under  different  object 
statistics.  The  maximum  Q algorithm  seems  to  compensate  optimally 
for  each  of  the  different  object  ensembles.  Again,  the  maximum  value 
of  Q is  approximately  one  unit  higher  than  tl^e  minimum  mean-squared- 
error  value. 

In  Fig.  8,  Q is  plotted  for  the  mean-only  nonstationarity.  Part 
(a)  has  A = 1,  and  (b)  has  A = 20.  The  blurred  images  in  each  part 
have  been  restored  using  the  maximum  Q filter  and  the  space-variant 
and  space-invariant  minimum  mean-squared-error  filters.  It  can  be  seen 
from  this  figure  that  there  is  very  little  difference  in  performance 
between  the  space-variant  and  space-invariant  minimum  mean-squared -error 
filters.  This  difference  does  not  depend  significantly  on  the  value  of 
the  object  center  mean-to-variance  ratio  A,  whereas,  in  earlier  figures 
of  space-variant  blurs,  the  loss  in  performance  suffered,  when  using  the 
space-invariant  filter,  depended  strongly  on  the  value  of  A.  As  before , 


the  maximum  value  of  Q is  relatively  easy  to  compute  and  parallels  the 
minimum  mean -squared -error  Q.  Something  not  seen  in  this  figure  is  the 
extreme  difficulty  in  computing  a space-invariant  maximum-Q-algor ithm 
filter  when  the  variance  Is  stationary. 

Figure  9 shows  Q plotted  for  images  restored  with  four  filters. 
These  are  the  space-variant  and  space-invariant  minimum  mean-squared- 
error  filters,  the  maximum  Q filter,  and  the  filter  derived  using  the 
maximum  Q algorithm  with  a space-invariant  filter.  Part  (a)  has  vari- 
ance-only nonstationary  and  part  (b)  has  mean  and  variance  nonstationary. 
Note  here  that  the  crosses,  which  are  the  results  of  space-invariant  fil- 
tering, show  quality  comparable  to  the  two  space-variant  filters.  This 
indicates  a computationally  workable  method  of  doing  useful  filtering  on 
nonstationary  blurred  image  ensembles.  Here,  the  maximum  Q filter  and 
space-invariant  maximum-Q-algor ithm  filter  are  both  easy  to  calculate 
and  both  parallel  their  respective  minimum  mean-squared-error  filters  in 
performance . 

A further  look  at  Fig.  9 reveals  two  related  points.  First,  as  might 
be  expected,  the  addition  of  the  nonstationary  mean  condition  of  Fig.  8 
to  the  nonstationary  variance  of  Fig.  9(a)  leads  to  a drop  in  both  the 
maximum  Q and  the  minimum  mean-squared-error  Q.  The  levels  of  Q for 
the  space-invariant  maximum-Q-algorithm  filter  and  the  minimum  mean- 
squared-error  filter  are  increased  by  allowing  t.he  mean  to  become  non- 
stationary. This  is  probably  due  to  "cancelling"  errors,  another  acci- 
dental improvement. 
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(d ) Filter  Error  Analysis 


Often,  comments  are  made  to  the  effect  that  there  are  many  theoret- 
ically good  ways  to  restore  degi*aded  images,  but  they  often  do  not  work 
well  due  to  one's  inability  to  estimate  the  object  statistics  and  blur 
impulse  responses  accurately.  We  have  therefore  done  the  following 
short  study  concerning  the  performance  of  our  filters  when  they  are 
derived  from  incorrect  statistical  or  blur  response  assumptions. 

The  curves  of  Fig.  10  show  Q for  a space-variant  blur  with  wc  = 8 

and  W =16.  The  object  ensemble  is  stationary.  Part  (a)  shows  the 
e 

performance  of  a space-variant  minimum  mean-squared-error  filter,  and 
part  (b)  is  that  for  a maximum  Q algorithm  filter.  The  three  curves 
in  each  part  include  restoration  with  a properly  constructed  filter  and 
two  restorations  with  filters  constructed  with  ±10$  errors  in  both  W 

c 

and  W . 
e 

In  Figs.  11  and  12,  a space-invariant  blur,  = 8,  is  used  with 
an  object  nonstationary  in  both  mean  and  variance.  In  both  figures,  the 
filters  used  are,  as  in  Fig.  10,  one  accurate  filter  and  one  filter  with 
approximately  10$  errors  in  some  parameter.  In  Fig.  11,  the  error  fil- 
ters have  a ±9$  error  in  the  object  width  W . In  Fig.  12,  there  is  a 

o 

±10$  error  in  the  blur  width  W^. 

As  can  be  seen  from  these  figures,  an  approximately  10$  error  in 
the  estimated  object  statistics  produces  only  a few  percent  reduction  in 
image  quality,  whereas  roughly  the  same  error  in  estimated  blur  produces, 
at  least  for  large  N, • a dramatic  loss  of  quality.  The  basic  reason  for 
this  dramatic  difference  is  as  follows.  An  error  in  object  statistics 
affects  only  the  point  where  the  "cutoff"  of  the  Wiener  filter  occurs. 
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Fig.  10.  Q for  stationary  object  and  space-variant  blur, 
Wc  = 8 and  We  =16,  for  mismatched  restoration  filter. 
Part  (a)  minimum  mean-squared-error  space-variant  filter 
and  part  (b)  maximum  Q filter.  In  both  parts,  filter 
parameters  are:  (1)  \VC  = 8 and  '.Ve  = 16;  (2)  Wc  = 8.8 
and  W = 17.6;  (3)  W = 7.2  and  W =14.4. 
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Fig.  11.  Q for  space-invariant  Wb  = 8 blur  and  nonsta- 
tionary object  mean  and  covariance,  WQ  = 34,  for  filter 
mismatched  in  object  width.  Part  (a)  space-variant  mini- 
mum mean-squared-error  filter,  and  part  (b)  maximum  Q 
filter.  In  both  parts,  filter  parameters  are:  (1) 

W =34;  (2)  W = 37.5;  (3)  W = 30.5. 

o o o 
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Fig.  12.  q for  space-variant  Wb  = 8 blur  and  nonstation- 
ary mean  and  covariance,  WQ  =34,  for  filter  mismatched 
in  blur  width.  In  part  (a)  minimum  mean -squared -error 

space-variant  filters  ar.d  part  (b)  maximum  Q filters. 

In  both  parts:  (1)  W =8;  (2)  W =8.8;  (3)  W =7.2. 

nb  b 


This  effect  should  be  most  pronounced  for  small  N where  the  signal-to- 


noise  ratio  is  small.  The  only  error  this  produces  is  either  a loss  of 
signal  or  increase  in  noise  for  those  spatial  frequencies  near  the  cut- 
off. 

An  error  in  estimating  the  blur,  on  the  other  hand,  means  that, 
when  N becomes  large  and  the  Wiener  filter  approaches  an  inverse  fil- 
ter, one  is  only  removing  a part  of  the  existing  blur,  or  possibly  adding 
additional  bur.  This  effect  will  be  known  to  anyone  familiar  with  in- 
verse filtering  techniques. 

This  result  is  somewhat  encouraging,  however,  since,  in  the  atmos- 
pheric and  compensated  imaging  systems  of  interest  here , the  amount  of 
blur  present  can  be  measured  with  reasonable  accuracy  while  the  object 
statistics  are  relatively  unknown  a priori. 
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APPENDIX  A 


In  this  appendix,  we  derive  a set  of  conditions  which  a linear 
filter  must  satisfy  to  insure  at  least  a locally  maximum  value  of  image 
quality  Q.  We  begin  with  the  definition  of  Q in  the  discrete  formu- 
lation 


where 


and 


i = N2  Tr I 


[H  ] [B  ] ] [Bj 

o 


,t 


[H] 
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(A  .1) 


(A. 2) 


[B]  o 


-2  /»  v /v  t \ t 

+ N [B][St  ][B]  ChF 
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2kN2[Sl[<l  l[B]tCHlt  + k2N2  [S]  [3t  HS]* 
o o 


(A  .3) 


The  procedure  we  shall  follow  here  is  to  assume  a form  for  an  ar- 
bitrary filter  [H]  which  includes  a scalar  parameter  t.  We  then  find 
the  first  and  second  partial  derivatives  of  Q with  respect  to  t and 
Jose  appropriate  conditions  on  these  to  guarantee  a local  maximum  in 
Q.  In  order  to  simplify  the  maximization  procedure,  the  quantity  log  Q 
will  be  maximized  instead  of  Q itself.  This  is  permissible  for  Q ^ 0 
since  the  logarithm  is  a monotonically  increasing  function. 

If  [H,,,.]  is  the  restoration  filter  which  produces  a maximum  image 
MQ 

quality  Q , then  any  arbitrary  filter  [h]  may  be  written 
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[H]  = Dl,„J  + t [g] 

MQ 


(A. 4) 


where  t is  a scalar  and  [G]  is  some  matrix.  A necessary  condition  for 
an  extremum  or  point-of-inf lection  in  Q is 


d log  Q 

5t 


= 0 vfc] 


(A  .5) 


t=0 


where  V means  for  all,  and  a necessary  condition  for  a maximum  is 


d log  Q 


dt" 


< 0 V [G] 


(A. 6) 


t=0 


Beginning  with  condition  (A. 5),  it  is  first  useful  to  note  that 

log  Q = log~^  - log  6 . (A. 7) 


Substituting  (A. 4)  into  (A. 3)  and  (A. 2),  and  taking  the  derivative  of 
the  logarithms  gives 


d log  Q (t ) 

5t 


2N2  Tr  (DHmq]  [B]  [SH  [B]1  [G]t  + t [G]  [B]  [Kq]  [B]t  tG]1) 


rf(t ) 


2Tr  ([H  1 [D]  [G]t  + t[G][D][G]t  - k>'2[G][B][«  ] [S]*) 

\ MQ O / 


6(t ) 


(A. 8) 


where  Q(’),  4(0,  and  6(0  indicate  the  values  of  Q,  <5,  and  6, 
respectively,  with  t = ••  Letting  t approach  zero  gives 
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j 


d log  Q (t ) 


dt 


2N2  Tr  (CH^l  [B 1 RU  CB]t  [G]*) 


t=0 


HQ 

^(t  = 0) 


2Tr([H  HdHG]  - kN2[G][B][»o][S]t) 


S(t  = 0) 


= 0 


V[G] 

(A. 9) 


This  resulting  condition  (A. 9)  will  be  satisfied  if 


°WI 


' /V  - 

[B]  o 


/Q  - 1' 

* l , 


CbHajCB]*  I - = CO] 


(A. 10) 


where  [0]  is  a matrix  with  all  zero  elements  and  again 


Q = q (0 ) = 

max  ^ J 5(0) 


(A. 11) 


This  analysis  eads  to  the  following  expression  for  the  maximum  Q fil- 
ter: 


LH  ] = kNfsm  JCB]* 
MQ  o 


CB1CK  ](£]* 

o 


(A. 12) 


Here,  it  must  be  remembered  that  in  most  cases  is  not  known  with- 
out first  knowing  • This  filter  is  sure  to  produce  a local  extre- 
mum or  point  of  inflection  in  Q. 

In  order  to  insure  a local  maximum,  we  need  condition  (A. 6).  Taking 
the  derivative  of  (A. 8)  at  t = 0 gives 


36 


d2  log  Q(t) 

N2  Tr(2[G][B][^ol(B]t[G]t) 

N4  Tr  (2  [H,„J  [B]  [S  ][B]^’tG]t) 

\ MQ  0 / 

>,  2 

ot 

<J(0) 

t=0 

<J2(0) 

Tr  ^2  [G  ] [D  ] [G  ] 1 ) Tr(2[H^1Q 

] [D] [G]1  - 2kS2lG)  [B]  [aQ] rs]t) 

6(0)  + 

62(0) 

(A. 13) 


After  some  manipulation,  we  have 


2NQmax  Tr 


o log  Q 


dt 


[B]  o 


ax 


-1\ 


max 


StBltS  ItBl'llOl* 


\ 

t 

/ 


6(0) 


t=0 


/ Matrix  \ 
lexpressiony 


l~gg  Q 
dt 


(A. 14) 


where  we  have  recognized  the  last  factor  as  the  first  derivative  at  t-0. 
Finally,  using  (A. 6)  with  (A.  14)  gives  the  condition 


2NQmax  Tr  [G] 


• /\  ^ 
[B]  o 


^Q’iax  Wb][«o][b]1:  ](G]' 


ax 


6(0) 


< 0 V[G]  . 


(A. 15) 


This  condition  is  insured  if  and  only  if 


t Qmax  ~ 1 5 [ill  Ub)1  > 


(0] 


(A. 16) 
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Thus,  if  conditions  (A. 9)  and  (A. 15)  are  met,  we  can  be  assured  of 
a filter  which  will  produce  a local  maximum  in  Q. 
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